%% Segmentacja obrazow za pomoca algorytmu HMM
% Katarzyna Rusinowska
% Marcin Krzewski
% Jakub Sitnicki

%% Reading images
clear all;
close all;

% wielkosc blokow
A = 32;
B = 32;

% wczytanie obrazkow
limg(1, :, :) = imread('D1.gif');
limg(2, :, :) = imread('D3.gif');
limg(3, :, :) = imread('D6.gif');

[I, M, N] = size(limg);

% struktura przechowujaca obliczone prawdopodobienstwa
probi(1:I) = struct('Pd_q', zeros(32, 1), 'pm_q', zeros(2, 256));

%% Learning HMM
for i = 1 : I
    img(1:M, 1:N) = limg(i, 1:M, 1:N);

    for j = 1 : 4
        x = mod(ceil(rand(1) * M), M/64) * 64 + 1;
        y = mod(ceil(rand(1) * N), N/64) * 64 + 1;

        crp_img = imcrop(img, [x, y, A-1, B-1]);

        Lambda = hmm(crp_img);
        
        Pd_q_tmp = stp(Lambda, true);
        pm_q_tmp = op(Lambda, crp_img, true);

        probi(i).Pd_q = probi(i).Pd_q + Pd_q_tmp;
        probi(i).pm_q = probi(i).pm_q + pm_q_tmp;
%         figure;
%         imshow(crp_img);
    end
    probi(i).Pd_q = probi(i).Pd_q ./ j;
    probi(i).pm_q = probi(i).pm_q ./ j;
end

%% Classification
seg_img = imread('img1.gif');
[M, N] = size(seg_img);

img = zeros(M, N);

 for m = 1 : M/A
     for n = 1 : N/B
        seg = imcrop(seg_img, [(n-1).*B+1, (m-1).*A+1, B-1, A-1]);
        
        Lambda_seg = lambda(seg); % hmm(seg);

        [Lc1, Ls1] = pl(Lambda_seg, seg, probi(1).Pd_q, probi(1).pm_q);
        [Lc2, Ls2] = pl(Lambda_seg, seg, probi(2).Pd_q, probi(2).pm_q);
        [Lc3, Ls3] = pl(Lambda_seg, seg, probi(3).Pd_q, probi(3).pm_q);

        % porownanie L1 < L2
        if comp(Ls1, Lc1, Ls2, Lc2) & comp(Ls2, Lc2, Ls3, Lc3)
            img([(m-1).*A+1 : m.*A], [(n-1).*B+1 : n.*B]) = 250;
        elseif comp(Ls1, Lc1, Ls2, Lc2)
            img([(m-1).*A+1 : m.*A], [(n-1).*B+1 : n.*B]) = 150;
        else
            img([(m-1).*A+1 : m.*A], [(n-1).*B+1 : n.*B]) = 50;
        end
     end
 end

imwrite(img, 'seg2.gif');
